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Abstract 

One of the major challenges of contemporary mathematics is numerical solving of various problems for functional 
differential equations (FDE), in particular Cauchy problem for delayed and neutral differential equations. Recently 
large variety of methods to handle this task appeared. However, many recent papers contain ambiguities or mistakes 
in either posing Cauchy problem for FDEs or in method of solving it. To rectify this situation, we propose a correct 
and efficient semi-analytical approach for FDEs based on Taylor method in this paper. The idea is to combine the 
method of steps and differential transformation method (DTM). In the latter, formulas for proportional arguments and 
nonlinear terms are used. Validity and efficiency of the presented algorithm is demonstrated by examples on wide 
classes of FDEs with multiple constant, time-dependent and/or proportional delays including FDE of neutral type. 
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1. Introduction 

Many engineering and chemical processes as well as economic and biological systems are modeled by a set of 
nonlinear delay differential equations (FDEs). As examples, we mention models describing machine tool vibrations 
01, behaviour of the central nervous system in a learning process, species populations struggling for a common food, 
dynamics of an autogenerator with delay and second-order filter, evolution of population of one or more species etc. 
Further models and details can be found for instance in Hi or lINl . 

To study the local dynamics of FDEs, the nonlinearities are expanded in a Taylor series about a fixed point or a 
periodic solution. The latter case gives rise to linear and nonlinear FDEs, depending on the degree of the series, with 
time-periodic coefficients. Recent example of Taylor series approach applied to study the local dynamics can be found 
in 0]. 

In the last two decades, various methods such as homotopy analysis method (HAM), homotopy perturbation 
method (HPM), variational iteration method (VIM), Adomian decomposition method (ADM), and also variety of 
methods derived from Taylor series method such as Taylor polynomial method, Taylor collocation method and dif¬ 
ferential transformation method (DTM) have been considered to approximate solutions of certain classes of FDEs 
in a series form. However, in several papers, for example Evans and Raslan 01, initial problems are not properly 
defined. The authors used only initial conditions in certain points, not the initial function on the whole interval, thus 
the way to obtain solutions of illustrative examples is not correct. Moreover, formulas used in the calculations are very 
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complicated. Similar situation occured in paper Blanco-Cocom, Estrella and Avila-Vales @], where the algorithms 
are not able to approximate a solution starting with another initial function. Both mentioned papers are using ADM. 
Comparison of ADM and DTM, which we decided to employ, is done in J3]. 

The concept of differential transformation in the form we utilize was proposed by Zhou {8[| in 1986. He applied 
the method to solve linear and nonlinear initial value problems in electric circuit analysis. This method constructs 
a semi-analytical numerical technique that uses Taylor method for solving of differential equations in the form of a 
polynomial. It is different from high-order Taylor series method which requires symbolic computation of necessary 
derivatives of the data functions. This is not the only way how to employ Taylor approach to solving FDEs. One 
pos sibility is to start with characterictic quasipolynomials and proceed to polynomial quasisolutions as in jsj] and 
[ hJ. It is also possible to derive the formulas using more general approach involving Hilbert spaces, see El or El- 
The crucial idea of our Taylor series concept is to combine general method of steps suitable for Cauchy problems 
for FDEs with a Taylor series technique known as DTM. For more details on method of steps see e.g. monographs 
Kolmanovskii and Myshkis H, Hale and Verduyn Lunel |d or Bellen and Zennaro El- Our approach enables 
us to replace the terms involving delay with initial function and its derivatives. Consequently, the original Cauchy 
problem for delayed or neutral differential equation is reduced to Cauchy problem for ordinary differential equation. 
Also the ambiguities mentioned above are removed. Moreover, while ADM, HAM, HPM and VIM require initial 
approximation guess and symbolic computation of necessary derivatives and, in general, 72 -dimensional integrals in 
iterative schemes, presented algorithm is different: Cauchy problem for FDE is reduced to a system of recurrent 
algebraic relations. 

2. Preliminaries 

2.1. Problem Statement 

For the purpose of clarity, we consider the following system of p functional differential equations of 72 -th order 
with multiple delays ari(f),..., a, -(f): 

u ( '°(0 = f(t, u(0, u'(f), - - -, u (n-1) (0, ui(ari(0), u 2 (ar 2 (0), - - -, u r (a r (t))), (1) 

where u {n \t) = (uf\t),..., u^{f)) T = uf\t),...,u p \t)), k = 0,1 ,...,72 - 1 and f = {fi,...,f p ) T are p- 

dimensional vector functions, u,-(a,-(f)) = (u(a,-(f)), u'(a,-(f)), ■ - ■. are 0«; • /?)-dimensional vector functions, 

nii < n, i = 1,2,..., r, r e N and fj : [0, 00 ) x R" p x R°' p are continuous real functions for j = 1,2,...,/?, where 

r 

LO — Yj m i- 

i=l 

We consider three types of delays a,: 

1. a,(f) = t - r,-, where r, > 0 is a real constant (constant delay) 

2. a,(t) = q,t, where q ,- 6 (0,1) (proportional delay) 

3. a,(t) = t - Ti(t), where r,-(f) > 0 is a real function (time-dependent or time-varying delay) 

Let f - min inf(a;(t))|, m - max{mi,m 2 ,... ,m r ], hence t* < 0 and m < n. In case m = n system Q is of neutral 

\<i<r [ t>0 ’ 

type, otherwise it is a delayed differential system. 

If f < 0, initial vector function <t>(f) = (<^i(f),..., cf> p (t)) T needs to be assigned to system (0} on the interval [t*, 0]. 
Further, for the sake of simplicity, we assume that <f>j{t ) e C n ([t*, 0]) for j = 1,...,/?. 

Consider system <|T]> subject to initial conditions 

u(0) = vo,u'(0) = vi,..., u (n_1) (0) = v„_i (2) 

and subject to initial vector function <t>(r) on interval [f, 0] such that 

0(0) = u(0),..., O (n_1) (0) = u (n_1) (0). (3) 

We consider Cauchy problem <[T]), Q and ([3]) under the following hypotheses: 

(HI) Suppose that Cauchy problem ([T]). (0 and 0 has a unique solution on some interval [0, T*]. 
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(H2) If aft ) = q,t and »;,■ = n in fj for some i G {1,..., r] and j G {1that is, y'-th equation is neutral with 
respect to proportional delay, we assume that uffaft)) = 0 for Z G {1,... ,p}, l + j. 

(H3) Suppose that the functions fj, j - 1 are analytic in [0, T *] x R np x K° p . 


Remark 1. Hypothesis (HI) is guaranteed for instance if the functions fj are continuous with respect to t on [0, 7”"] 
and Lipschitz continuous with respect to the rest of the variables on R" p X R up , the functions fj, tp '-,..., are Lip- 
schitz continuous on [t* , 0] and a ,• are Lipschitz continuous on [0, T*\ For more details, we recommend publications 

MMorM- 

Remark 2. Hypothesis (H2) is included since if it does not hold, the existence of unique solution of Cauchy problem 
could be violated. 


2.2. Differential Transformation 

Differential transformation of the k-th derivative of function u(t) is defined as follows: 


m) - h 


d k u(t) 

dt k 


t=t 0 


(4) 


where u(t) is the original function and IJ(k) is the transformed function. Inverse differential transformation of U(k) is 
defined as follows: 

oo 

u(t) = Y J u m- t of, ( 5 ) 

it=0 


hence 


OO 


u(t) = 'Yj 

k=0 


(f ~ h)) k 
k\ 


d k u(t) 

dt k 


t=t 0 


( 6 ) 


In fact, formula ([6} indicates that the concept of differential transformation is derived from Taylor series expansion. 
However, DTM is not evaluating the symbolic derivatives: instead of this, relative derivatives are calculated by an 
iterative way which is described by the transformed form of the original equation. 

The following well-known formulas for to = 0 can be easily derived from definitions ©. (0: 


Lemma 1. Assume that F(k), H(k) and U(k) are differential transformations of functions f(t), h(t ) and u(t), respec¬ 
tively. Then: 


if m 
lf/(0 
If/W 

if/(0 


t n , then F(k) = 6(k - n), where his the Kronecker delta. 


e At , thenT’(k) = —. 

k\ 

d n u(t) (k + n)\ 

-^,then F(k)= [ - 1 ^U(k + n). 

k 

u(t)h(t), then F(k) = ^ U(l)H(k - /). 

1=0 


Similarly transformation formulas for proportionally delayed arguments can be proved, as in & 
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Lemma 2. Assume that W(k),U(k),Ui(k) are the differential transformations of the functions w(t),u(t),Uj(t ) and 
q, qt e (0,1), i = 1,2 then: 


If w(t) = u(qt), then W(k) = q k U(k). 

k 

If wit) = u\{q\f)U2{q 2 t), then W(k) = ^ q\q k ff l Ufl)U 2 {k - /). 

1=0 


d'"u(qt) (k + m)\ k 

If w(t) = ——— .then W(k ) =-—— q k U(k + m). 


d{qt) n 


k\ 


The main disadvantage of standard form of DTM is that it cannot be used directly for equations with nonlinear 
terms containing unknown function u(t). For example, for fin) = V1 + if or fin) — e sin “ etc. 

However, differential transformation of components containing nonlinear terms can be calculated using the so- 
called Adomian polynomials A„ in which each solution component u, is replaced by the corresponding differential 
transformations component U(i), i = 0,1,2,... (see Q). If F(k) is the differential transformation of a nonlinear 
term f(u). then 


F(k) = A„(f/(0), 1/(1),. ■ U{n))6(k - n) = A k (U{ 0), 17(1),. ..,U(k)) 

n =0 


1 d k 
ki.dtk 




\i =0 


t =0 


k> 0. 


For illustration, we give several components explicitly: 


7 ( 0 ) - num, 

7(i) - u(i)f(um, 

7(2) = U(2)f'(Um+^U 2 (l)f"(Um 

7(3) = f/(3)/'(f/(0)) + U(l)U(2)f"(Um + i(/ 3 (l)/"'([/(0)). 


(7) 


3. Main Results 

Recall system (Q} 

u M (t) = f (r, u (t), u'(t), u ( ' !_1) (0, ui(ai(0), u 2 (a 2 ( 0 ), • • ■, u r (a r (t))), {Q> 

with initial conditions 

u(0) = v 0 , u'(0) = Vi,...,u ( ” _1) (0) = v„_i © 

and initial vector function <F(f) on interval [f, 0] satisfying 

4>(0) = u(0),..., ® (n_1) (0) = u ( " _1) (0). © 

First we apply the method of steps. We substitute the initial function <F(f) and its derivatives in all places where 
unknown functions with constant or time-dependent delays and derivatives of that functions appear. Then system Q 
changes to a system of ordinary differential equations or differential equations with proportional delay, depending on 
whether system © contains terms with proportional delays or not. For example, if aft) = q\t, a 2 {t) = t — T 2 (t) and 
a ft) - t - Tj, aft) = t — T4 are constant delays, applying the method of steps transforms system © to system 

U (n) (0 = fit, u(0, u'(f), . . . , U (n ~ l \t), m ( 9 i0, ®2(f - T 2 (0), (t - Ts), ®4(f - r 4 )), (8) 
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where mfoiO = (u(q 1 t),u'(q 1 t),...,u^\q l t)), ® 2 (f - r 2 (0) = Mt - r 2 (0), ■ ■ ■, ® (m2 \t - T 2 (f))), ®/(* - r,-) = 
(Off - Ti), O'ff - r,-),. .., - t,)), i = 3,4, m, < n, l = 1,2,3,4. 


Next we transform initial conditions (0. Following definition of differential transformation, we derive 

U(*) = ^u w (0). 
k\ 

Applying DTM to system ([8]) we get system of recurrent algebraic equations 

U (k + n) = r(k, U(0), U(l),..., U (k + n- 1)). (9) 

Using transformed initial conditions and then inverse transformation rule, we obtain approximate solution of system 
0 in the form of infinite Taylor series 

oo 

(t=0 

If t* = 0, which means that all delayed arguments are proportional, the approximate solution u(f) is valid on the 
intersection of its convergence interval with [0, T*]. Otherwise, if f < 0, we denote t aj = inf jf : aft) > 0} and 
t„ = min {t a . : t n 4 0). Then the approximate solution u(f) is valid on the intersection of its convergence interval and 

l</<r ' ' 

the interval [0, t a \ n [0, T*], and u(f) = Off) on the interval [f*, 0]. 


Now we are prepared to formulate and prove the following main result: 

Theorem 1. Let the hypotheses (HI), (H2) and (H3) hold. Then the unique solution of Cauchy problem consisting of 
system 0 together with conditions 0 and © eventually, is locally approximate by using combination of method of 
steps and differential transformation method or differential transformation only. 

Proof. If all the assumptions are fulfilled, then either we are dealing with the system of equations with proportional 
delays only, or we can pass to such system or system of ODEs using the method of steps. Validity of this approach 
can be verified for instance in 0, J1 or 0. 

The obtained system has a unique solution subject to initial conditions 0, either on [0, T*] or on [0, t a ] n [0, T*]. 
Then, since the righthand side of the system is analytic, differential transformation leads to solving a system of 
difference equations instead of solving the system of differential equations. 

Using the transformed initial conditions, a finite part of the sequence which is a solution of difference equation 
is calculated. From the definition of differential transform we see that the finite sequence we calculated represents 
coefficients of Taylor polynomial of the solution of Cauchy problem. According to the Theorem 2.2.1 in 1 1311 . this 
polynomial is a local approximation of the unique solution which is valid on some interval [0, b], 0 < 6 < T* where 
T* = min{f«, T*}. 

Corollary 1. For the unique solution of Cauchy problem 0, 0 and © obtained by Theorem\J\ we have the following 
error estimate on [0, b]: 

K 


where N is the order of Taylor polynomial approximation and K is a bound for (N + 1 )-th derivative of the solution 
on [0, b]. 

Proof. Since the functions fj are analytical in [0, T*] x R np x R ojp , without the loss of generality we can suppose 
that using Theorem [T| we are able to approximate the unique solution by Taylor polynomial 

N 

u N (t ) - U(k) f* 

k =0 


5 



in [0,5]. The remainder of Taylor series is given by 


it — un —Rn — 


1 


j^+i 


u(t) 


(N + 1)! dt N+1 


t=h 


for some t\ such that 0 < t\ < t. Since (N + l)-th derivative of u(t) is bounded for t 6 (0, <5], i.e. 


d N+1 u(t ) 
dt N+l 


< K 


for some nonnegative constant K, then the maximum error for u N (t) in this interval can be estimated from this remain¬ 
der term as 


K 

6N ~ (ivTTj!' 


4. Applications 

Example 1. Consider the nonlinear system of delayed differential equations 

t _ 2 

hi | — U 2 5 



M 3 = e 2 'ii 2 (t) + 9e 2, u 3 



subject to the initial conditions 


M] (0) = 1, M2(0) = 1, m 3 (0) = 0. 


Applying DTM to system (ITOt we get recurrent system 

k 

(k+V)U l (k+\) = Y J U2{l)U2{k-D, 

1=0 


(k + \)JJ 2 (k + 1) = 
(k+l)U 3 (k+l) = 



U 3 (k-l). 


From initial conditions we have U i(0) = 1, U 2 (0) - 1, U 3 (0) = 0. Solving recurrent system (fl2T> we get 


k = 0 : t/i(l) = t/|( 0 ) = 1 , 

u 2 ( 1 ) = l -um =1 

f/ 3 (l) = f/ 2 (0) + 9t/ 3 (0) = 1. 


( 10 ) 


( 11 ) 


( 12 ) 


k= 1 : f/j(2) = i (t/ 2 (0)f/ 2 (l) + U 2 (l)U 2 (0)) = i 
U 2 ( 2) = |c/r(l) = ^ 

U 3 ( 2) = |(C/ 2 (1) + |tf 2 (0)) + ^ (3f/ 3 (l) + 2t/ 3 (0)) = 3. 
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Consequently for k > 3, we find 


U i(3) = i 




^2(3) - 4g , £/2(4) - 384 , U 2 (5) - 3840 ,-- 

rr /a\ _ 9 „ / / i\ _ 27 7 , /c 3 _ 81 


1/3(3) = §, 


t/ 3 (4) - f, 




Using inverse differential transformation we obtain a solution of ( ITOt . dTTb in the form 


„ , 1 , 1 3 1 4 1 4 

ui(t)-\ + t + —r + —r + —t + —t + 


2 ! 

1 1 


4! 


5! 


W2 ( 0 - 1 + -f + -f 2 + _r 3 + 


w 3 (7) — 7 + 3f 2 + 


2 9 2 

u • -t + 


1 4 

384 f + ' 
81 


27 3 01 4 

—7 + —7 + 
6 24 


(jLiiLliLC 

L 2! 3! 3! 4! 

\ 2! 3! 4! ) 


•= 1 + 2 + 


= 7 


If k —> 00 , then the series solution converges to the series expansion of the closed form solution 

«i(7) = e', u 2 (0 = e' l2 , m 3 (7) = te 3 '. 

The same initial problem has been solved by Saeed and Rahman El using Adomian decomposition method but 
authors obtained only approximate solution. 

Example 2. Consider the nonlinear system of delayed differential equations 

u'l = 2u\ j u 2 (t - 1) + u\ (?) - t 4 + It 3 - r - 4? + 4, 

«2 = U2 (^) Ul {^) + u ~ 2 ) “ Y 8 ~ 2t + 6 ( 13 ) 

with initial functions 


MO = 2 ? 2 , 

MO = t 2 (14) 


for 7 e [-2,0] and initial conditions 


u i(0) = (0) = 0, u\ (0) = 0j(O) = 0, 

U 2 (0) = 02(0) = 0, «2(0) = 0i(O) - 0. 

Using the method of steps we get 

u" = 2(7 - 1) 2 mi + u\(t) - t 4 + 2t 3 - f 2 - At + 4, 

Applying DTM to system {[6} we get recurrent system 

(k + 1 )(k + 2)U\(k + 2) = 2 ^ (6(1 - 2) - 26(1 - 1) + 6(1)) - Ui(k - l) + (k + \)U x (k + \) 

to l 2 / 

- 6(k - 4) + 26(k - 3) - 6(k - 2) - 46(k - 1) + A6(k) 


(15) 


(16) 


(17) 
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and 


k 

(k+ \)(k + 2)U 2 (k + 2) = ^ 
1=0 



U 2 (l)Ui(k - /) - —6(k - 4) + 26{k). 
to 


From initial conditions we have t/ 1 (0) = 0, (7i(l) = 0, 

U 2 ( 0 ) = 0, t/ 2 ( 1) = 0. Solving recurrent system (fl7T> . (fTSl i we get 


k = 0 : 2U\(2) = 2U\(0) + t/i(l) + 4 => U i(2) = 2, 
2(7 2 (2) = U 2 (0)Um + 2^ U 2 ( 2) = 1. 


( 18 ) 


A: = 1 : U i(3) - -Wl) + 2(7,(2) - 4(7,(0) - 4) = 0, 
o 

t/ 2 (3) = V(7 2 (0)(7,(0) + i(7 2 (l)(7 1 (0)) = 0. 

O 2 3 

k - 2: (7,(4) = I |(^(7i(2) - £/j(l) + (7,(0)) + ^(7,(3) - ^ = 0, 

W) = i ^(7 2 (0)(7,(2) + 1(7 2 (1)(7 1 (1) + i(7 2 (2)(7j(0)J = 0. 

Now it is obvious that for k > 3 all coefficients £7 1 (A" + 2), U 2 (k + 2) are zero. 

Using inverse differential transformation we get 

it 1 (0 = 2f 2 , w 2 (f) = f 2 

which is exact solution of initial problem (fl3l i. (fTTb . (fl5l> . 


Example 3. Consider the nonlinear system of neutral delayed differential equations 

w'i" = u'"(t - 2)u x ^ j + ^m 2 + u' 2 (t - 

u 2 ~ u 2 + l, 2^ ~ !) w i (t) ( 19 ) 

with initial functions 

M0 = e', 

M0 = t 2 ( 20 ) 

for t e [-2,0] and initial conditions 

u 1 (0) = 1, m'j(O) = 1, m"(0) = 1, 

u 2 ( 0) = 0, u 2 (0) = 0, m"(0) = 2. (21) 


Using the method of steps we get 

m'j" = e (f_2) Mi j + +yjui + 2t — e~', 

«'"=«'"(^) + 2(f- 1 ) M ,(^). (22) 

System (l22l) can not be solved using classic DTM with respect to nonlinear term f(u) = JW[, therefore we will apply 
modified Adomian formula for differential transformation components to nonlinear term f(u). At first, applying DTM 
to system (l22l> we get recurrent system 


(k + 1 )(k + 2){k + 3)U\(k + 3) = e~ 2 



U\(k-l) + F\(k) + 26(k) - 


(-!/• 

kl 


(23) 
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where F\(k) is transformed function of f(u) = According to formula ([7]) and from transformation of initial 

conditions 


um = i, t/id) = i, u 1(2) = 1 

U 2 ( 0) = 0, U 2 ( 1) = o, U 2 ( 2) = 1, 


we get 


^i(O) 

Fi(l) 

F i(2) 




f/f(0) = 1, 

2 £/i(l) _ 2 

3 WM ~ 3’ 

2 £/i(2) _ 1 

3 WT(0) ” 9 


£/?(!) 


5 

9’ 


Solving recurrent system (l23l >. (l24l > we have 

k = 0 : C/r(3) - + F,(0) + l) = 

t/ 2 (3) = ^(-2t/ 1 (0)) = -i 

k= 1 : C/i(4) = + t/1(0)) + 1) + 1) = 

t/ 2 (4)=^(2t/ 1 (0)-^t/ 1 (l)) = i 


fc = 2: 


c,,<5) = ^(5' /l<2) + 5 t/ ' <1) + ^'(O)) + s( f ' ,<2) " ’) 

1/2 2 v 4 

C/ 2 ( 5 ) = —(-E/i(l) - — t/ 1 ( 2 )) =-. 

45^3 9 ' 405 


17e" 2 - 8 
1080 


Using inverse differential transformation we get for initial problem (IT9t . ( 12Oi l. (I2H only approximate solution 


— 1 + t + t + 


2 + e- 


-t + 
4 


4e- 2 + 5 


72 


-t + 


\le~ 2 - 8 
1080 


-f + ..., 


u 2 (t) — It 1 — — T + — d 4 -- 1 + ... 

w 3 9 405 


5. Conclusion 

• We conclude that combination of method of steps and differential transformation method (DTM) presented in 
this paper is a powerful and efficient Taylor series technique suitable for numerical approximation of a solution 
of Cauchy problem for wide class of nonlinear functional differential equations, in particular delayed and neutral 
differential equations. No discretization, linearization or perturbation is required. 















• There is no need for calculating multiple integrals or derivatives and less computational work is demanded 
compared to other popular methods (Adomian decomposition method, variational iteration method, homotopy 
perturbation method, homotopy analysis method). 

• Using presented approach, we are able not only to obtain approximate solution, but even there is a possibility 
to identify unique solution of Cauchy problem in closed form. 

• A specific advantage of this technique over any purely numerical method is that it offers a smooth, functional 
form of the solution over a time step. 

• Another advantage is that using our approach we avoided ambiguities, incorrect formulations and ill-posed 
problems that occured in recent papers. 

• Finally, a subject of further investigation is to develop the presented technique for system ([]} with state depen¬ 
dent delays and for partial differential equations. 
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